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BLOCK  20.  ABSTRACT  (cont) 

Over  150  simulation  runs  were  performed  over  such  terrain  having  added  Gauss- 
Markov  noise  (independent  of  the  terrain  process)  with  6  different  noise  levels. 
Results  of  the  simulation  are  presented,  giving  performance  as  a  function  of 
signal/noise  ratio  and  as  a  function  of  terrain  roughness.  Measured  miss- 
distances  are  also  presented.  The  simulation  also  includes  the  performance  of 
a  nearest-prototype  classifier  based  on  mean  absolute  distance  (the  first 
Minkowski  metric) ,  which  is  the  actual  metric  employed  in  TERCOM. 

The  following  conclusions  are  reached:  (1)  the  analysis  offers  a  useful  descrip¬ 
tion  of  actual  TERCOM  performance  on  both  artificial  and  real  terrain;  (2)  the 
current  model  is  different  from,  and  more  convenient  to  use  than  the  model 
previously  available;  (3)  the  TERCOM  classifier  should  be  slightly  modified 
with  respect  to  the  mean-removal  technique  in  order  to  significantly  improve 
the  theoretical  and  the  actual  performance;  (4)  TERCOM  is  an  effective  naviga¬ 
tion  system  provided  that  the  terrain  possesses  certain  statistical  charac¬ 
teristics;  and  (5)  the  necessary  statistical  characteristics  of  the  terrain 
can  be  summarized  in  two  parameter  ratios,  oz/o t  and  0z/°tl»  which  relate  to 
terrain  and  noise  variance,  and  which  are  defined  in  the  report. 
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SUMMARY 


PROBLEM 

An  analysis  of  the  Terrain  Contour  Matching  Navigation  System  (TERCOM) 
was  performed  by  the  system's  manufacturer#  LTV  E-Systems,  Enc.  An  inde¬ 
pendent  evaluation  of  that  analysis  and  of  TERCOM  w as  desired  for  several 
reasons?  (1)  some  doubt  existed  about  initial  assumptions;  (2)  that  analy¬ 
sis  was  performed  for  a  mean-square-distance  (MSD)  classifier  with  continu¬ 
ous  input,  but  a  mean-absolute-distance  (MAD)  classifier  with  discrete 
input  is  used  in  the  actual  hardware;  and  (3)  the  utility  of  the  single 
terrain  parameter,  cz,  suggested  by  E-Systems  as  being  predictive  of 
TERCOM  performance  seemed  limited  because  it  varies  widely  over  terrain 
areas  of  interest. 

APPROACH 

The  approach  was  twofold: 

<1)  Starting  with  basic  assumptions,  derive  an  error  model  for 
TERCOM.  Concentrate  on  the  MSD  ►classifier  because  of  its  mathematical 
tractability,  but  treat  the  MAD  classifier  also  as  far  as  possible. 

(2)  Determine  by  computer  simulation  the  performance  capabilities 
of  MSD  and  MAD,  using  both  artificial  Gaussian  terrain  and  real  terrain 
samples.  Compare  TERCOM  performance  over  artificial  statistically  control¬ 
led  terrain  with  the  performance  over  real  terrain  to  yield  the  parameters 
necessary  for  a  predictive  performance  model. 

RESULTS 

A  mathematical  model  for  TERCOM  "false- fix"  probability  was  derived 
that  differs  from  the  original  LTV  E-Systems  model  in  certain  assumptions. 
False  fix  probability  is  presented  as  a  finite  sum  having  parameters  that 
depend  on  sample  size,  terrain  and  noise  correlations,  and  the  signal-to- 
noise  ratio. 

The  computer  simulation  revealed  that  a  dramatic  increase  in  the 
probability  of  a  correct  fix  could  be  achieved  by  a  simple  change  in  the 
normalization  procedure  for  the  MSD  or  MAD  computations. 

It  was  found  that  the  simulated  performance  of  TERCOM  on  real  terrain 
was  the  same  as  on  artificial  terrain,  generated  from  a  Gaussian  process, 
when  a  parameter  ratio  a z/ot  was  the  same  for  both  terrains.  This  per¬ 
formance  was  specified  by  plots  of  probability  of  correct  fix  vs  the 
parameter  ratio  a z/0n*  A  family  of  performance  curves  was  generated  by 
varying  the  parameter  ratio  oz/Of 
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The  error  model  developed  herein  is  a  good  predictor  of  simulation 
results,  but  a  rather  poor  predictor  of  performance  on  actual  terrain. 

This  is  believed  to  be  principally  due  to  the  assumption  of  stationary 
statistics  for  the  terrain.  While  a  Gauss-Markov  uvdel  may  be  valid  for 
local  terrain  regions,  the  Gauss-Markov  statistics  vary  widely  over  real 
terrain,  even  over  physical  extents  represented  by  the  reference  arrays 
encountered  in  TERCOM.  Both  the  E-Systems  error  model  and  the  present 
one  offer  insight  into  TERCOM  performance,  and  both  would  undoubtedly  be 
improved  by  considering  the  nonstationarity  of  the  terrain  statistics. 

The  present  model  is  not  as  computationally  involved  as  the  E-Systems 
model,  however,  and  may  therefore  be  preferable. 

CONCLUSIONS 

TERCOM  performance  can  be  significantly  improved  by  means  of  a  minor 
modification  in  the  mean-removal  procedure. 

Present  error  models  of  TERCOM  suffer  from  an  assumption  of  stationary 
statistics.  This  assumption  is  not  justified  on  real  terrain  except  when 
precautions  are  taken  to  insure  its  validity. 

A  parametric  family  of  performance  curves  can  be  generated  by  simula¬ 
tion  of  TERCOM  flights  over  modified  Gaussian  terrain.  Once  such  a  family 
of  curves  is  generated  for  a  particular  terrain  sample  size,  all  performance 
characteristics  of  interest  sxe  determined.  It  still  remains  to  compare 
these  results  with  actual  flight  test  data. 

TERCOM  is  an  effective  system  provided  that  reference  arrays  can  be 
obtained  that  produce  acceptable  performance.  Such  reference  arrays  are 
characterized  herein,  but  their  availability  in  the  real  world  was  not 
examined. 
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PREFACE 


The  work  discussed  in  this  report  was  performed  under  Project  7233, 
"Applications  of  Biological  Principles  as  Solutions  to  Air  Force  Needs 
in  Signal  Processing  and  Information  Handling,"  Task  7233-05,  "Automatic 
Pattern  Recognition  for  Air  Force  Systems,"  Work  Unit  7233-05-17,  "Appli¬ 
cation  of  Visual  Pattern  Recognition  to  USAF  Problems."  The  analysis  and 
simulation  were  performed  in  the  Mathematics  and  Analysis  Branch,  biody¬ 
namics  and  Bionics  Division,  6570  Aerospace  Medical  Research  Laboratory, 
during  the  period  19  October  1972  to  15  November  1973.  The  work  was 
performed  at  the  request  of  the  Subsonic  Cruise  Armed  Decoy  (SCAD) 

System  Program  Office  (ASD/RW86) . 

The  authors  thank  Dr.  O.H.  Tallman,  III,  Col,  USAF,  and  Dr.  H.L. 
Oestreicher  for  useful  discussions  and  encouragement  during  this  effort. 
Thanks  are  also  due  Capt  J.  Gobien  who  suggested  use  of  the  Toeplitz 
approximation  in  the  mathematical  analysis,  and  Capt  D.  Brungart 
provided  the  computer  graphics. 
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INTRODUCTION 


The  terrain  contour  matching  navigation  system  (TERCOM)  (Ref  1)  has 
been  proposed  as  a  self-contained,  autonomous  system  for  updating  an 
inertial  guidance  system  at  selected  checkpoints  en  route.  The  system 
is  based  on  a  pattern  classifior  that  may  be  termed  a  nearest-proto type 
classifier,  where  "nearest"  is  defined  i.n  the  first-Minkowski-metric  sense. 
In  order  to  obtain  a  clearer  understanding  of  the  capabilities  and  the 
limitations  of  the  system.-  an  analysis  and  evaluation  of  TERCOM  was  per¬ 
formed  by  this  laboratory  at  the  request  of  the  monitoring  agency. 

The  report  of  that  analysis  and  evaluation  is  organized  as  follows. 

The  first  section  presents  a  brief  discussion  of  the  performance  of  an 
error  model  for  TERCOM  that  was  developed  by  TERCOM' s  manufacturer.  Some 
theoretical  results  are  presented  as  suggestions  for  improving  the  model. 

A  brief  comparison  of  nearest  prototype  classifiers  based  on  the  first 
Minkowski  metric  (termed  MAU  classifiers)  and  based  on  the  Eucliuean 
metric  (termed  MSD  classifiers)  is  also  included.  The  next  section  pre¬ 
sents  a  discussion  of  the  simulation  procedure  and  the  simulation  results. 
The  last  section  presents  the  conclusions  and  recommendations. 

THEORETICAL  CONSIDERATIONS 
PERFORMANCE  OF  THE  TERCOM  ERROR  MODEL 

Because  the  theoretical  aspects  of  pattern  recognition  in  general, 
and  TERCOM  system  performance  in  particular,  are  quite  mathematical,  there 
must  be  some  fairly  strong  motivation  for  the  typical  systems  designer  to 
consider  them  in  depth.  The  initial  motivation  for  development  of  an 
error  model  is,  of  course,  to  predict  systems  performance:  such  predic¬ 
tions  are  essential  to  the  systems  buyer  and  user.  But,  if  the  TERCOM 
manufacturer  has  provided  an  error  model,  why  do  we  need  more  mathematical , 
theoretical  considerations? 

The  answer  is  obtained  by  constructing  a  scatter  plot  of  the  pre¬ 
dictions  of  the  manufacturer's  error  model  against  actual  flight  test  data. 
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Such  a  plot  would  reveal  that  predictions  do  not  correspond  well  with  the 
measured  performance.  The  Pearson  correlation  coefficient  computed  for 
such  a  scatter  plot  (p  =  0.617)  is  interprerable  as  meaning  that  the  error 
model  accounts  for  only  about  38%  of  the  variance  in  the  measured  per¬ 
formance.  This  leaves  about  62%  of  the  variance  unaccounted  for.  When 
wo  consider  the  cost  of  a  system  like  TERCOM,  there  is  clearly  strong 
motivation  to  reconsider  a  mathematical  error  model. 

GENERAL  CONCEPTS 

Let  x  =  (X0'X2'  *  •  •  ,x<j-i)be  t^ie  true-position,  stored,  sampled  terrain 

contour.  Assume  x  is  accurate  to  some  "mapping"  error  that  is  negligible 

compared  to  the  observation  noise.  Let  r  =  (rQ,r^, . . . ,r^  ^  =  x  +  n, 

where  n  =  (n  ,n, ,...,n.  ,)  is  the  observation  noise.  Let  y(x)  = 
o  1  d-1 

(yo,y^, . . . ,yd_^)  be  the  stored  sampled  terrain  contour  at  a  geographical 
distance  x  from  the  true  position,  as  shown  in  Figure  1.  Finally,  let  m 
denote  the  match-distance  (e.g.,  the  decision  statistic)  between  the 
vectors  r  and  y(f) ,  and  let  mQ  denote  the  special  case  T  =  0. 


END  OF 

TRACK,  X, 
a 


X 

START  OF 
TRACK,  Xj 


FLIGHT  PATH 


Figure  1.  Nome »ciature  and  notation  for  vector  samples 
from  a  stored  reference  array. 
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First  of  all  we  indicate  the  procedure  for  finding  the  probability 
that  m  <  mo  for  some  fixed  x.  This  will  show  us  what  distributions  we 
need  to  find  in  order  to  arrive  at  numerical  predictions. 

We  let  pT(m,mo)  denote  the  joint  probability  density  function  of  m 
and  mp.  Recall  that  m  corresponds  to  some  t,  and  the  function  depends 
(in  general)  on  this  x:  hence  the  subscript.  Then  we  compute 


Prob  {m  <  mg (mo)  = 


0 


PT(m,mo)  dm 


When  we  allow  mo  to  range  over  all  permissible  values  we  arrive  at 


r“  0 

Prob  {m  <  mn}  =  p  (m,mo)  dm  dmn 

JO  Jo  T 


We  now  make  our  first  assumption.  We  assume  that  for  sufficiently 
large  x,  the  distributions  of  m  and  mo  are  independent.  In  the  case  of 
the  Gauss-Markov  model,  which  we  introduce  shortly,  this  corresponds  to 
r  >  Lt,  where  J»T  is  the  correlation  length  of  the  terrain.  This  is 
equivalent  to  the  following  assumption: 


Prob  {m  <  mo)  = 


.00 

p(mo) 

0 


dm  dmo 


Now  we  want  to  know  the  probability  that  mo  is  smaller  than  m  for 
every  permissible  value  of  x.  Suppose  that  p^im)  is  the  same  for  each 
of  the  permissible  values  of  x  (but,  recall,  x  >  L^)  ,  and  that  the  inte¬ 
gral  above  has  been  evaluated  to,  say,  6.  That  is, 


Prob  {m  <  mo)  &  (3 

Then  Prob  {mo  <m}  =  l-  6=  cx 

Note  that  supposing  PT(m)  does  not  change  with  x  imposes  a  second  assump¬ 
tion:  that  the  terrain  statistics  are  stationary. 
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To  proceed,  consider  the  rank  ordering  of  three  observations  that 

correspond  to  mg ,  ,  and  m2 ,  where  m  and  m2  are  match  values  for  two 

distinct  t  >  L  .  Six  orderings  are  possible;  suppose  the  probability  of 

the  1  order  is  p.. 

1 1 


m0 

mj 

m2 

Pi 

mg 

m2 

mi 

P2 

mj 

mQ 

m2 

P3 

mj 

mg 

P4 

m£ 

m0 

mi 

PS 

m2 

mi 

mg 

P6 

We  note  that  Prob  {mg  <  m^}  =  P3  +  P4  +  P6  =  ct,  and  Prob  {mg  <  m2]  = 
P4  +  P5  +  p6  =  a.  It  is  a  simple  matter  to  compute 

Prob  (mg  <  mj  mg  <  m2}  -  P4  +  P6  •  We  further  note  the  constraint 

f 

'u  p.  =  1.  These  facts  together  imply  that 
i=l  1 

Prob  f r_ g  <  mx  m0  <  m2}  =  2o  -  1  +  Pi  +  P2 

Now,  as  a  approaches  one,  it  is  clear  that  pi  and  P2  must  approach  zero. 

Hence,  for  "large"  a, 

Prob  {mg  <  mi  mg  <  m2}  ~  2a  -  1 

As  a  approaches  zero.  Pi  +  P2  dominate.  If  we  suppose  Pi  =  P2  =  1/2,  we 
obtain 

Prob  (mg  <  mi  f~\  mg  <  m2}  =  a/2 

The  somewhat  surprising  thing  is  that  this  sort  of  combinatorial  argument 
can  be  generalized  to  rank  orderings  of  many  observations:  the  approxi¬ 
mations  remain  the  same. 

The  approximations  just  developed  are  fine  for  the  high  and  low 
ranges  of  a.  But  what  about  a  =  1/2?  One  case  where  o  equals  1/2  arises 

when  p^  =  1/N  for  all  i,  and  N  is  the  number  of  possible  rank  orderings. 
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If  M  is  the  number  of  observations  being  ranked,  then  N  =  Ml  and  Prob  {niQ 
is  the  minimum  observation}  =  1/M.  We  see  that  in  this  particular  case, 
the  approximation  for  "large"  a  gives  the  best  result  when  M  is  large. 

Using  the  preceding  combinatorial  argument  as  motivation,  we  choose 
the  following  approximation. 

Prob  {mo  <  all  m} 

^  max  { 2cx-l ,  0} 

=  max  {1-28,  0} 

F  fm  0 

=  max  { 1-2  p(mo)  p  (m)  dm  dniQ,  0}  (1) 

JO  JO  T 

This  expression  approximates  the  probability  of  correctly  assigning  r  to 
a  position  that  is  within  a  distance  of  L^  of  x. 

At  this  point  we  see  what  distributions  are  of  principal  interest: 
p(mo)  and  p^_  (m)  .  We  must  consider  the  definition  of  the  classifier’s 
match  function  before  we  can  proceed  further. 

THE  GAUSS-MARKOV  STATISTICAL  MODEL 

The  statistics  of  the  terrain  process  and  the  noise  process  are 
assumed  to  be  Gauss-Markov,  and  independent  of  each  other.  Hence  the 
probability  density  function  (using  the  notation  of  Ref  5) 

p-(5)  ■=  (2i5fd/2|AT|  /  exp  {-  j  aA^o1} 

governs  the  terrain  process,  and  the  probability  density  function 

P- (5)  =  (2u;  d/2iAJ  exp  {- r- 5A_1aT} 
n  n  i  n 

governs  the  noise  process.  The  covariance  matrix  for  the  terrain,  A  , 

tli  T 

has  x-j  elements  given  by 

02  -|i-j|AT 

T 
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th 

The  covariance  matrix  for  the  noise,  A  ,  has  i-j  elements  of  the  same 

n  J 

form,  but  with  different  parameters, 

a2  c-|i-l|An 

n 

The  parameters  and  are  called  the  correlation  lengths  of  the 
terrain  and  the  noise,  respectively. 

The  joint  distribution  of  x  and  y  (x)  has  a  similar  form,  but  the 
covariance  matrix  is  slightly  more  complicated.  Define 


A 


z 


A  ' 
xy  i 


where  Am  is  the  terrain  covariance  matrix,  as  before,  and  A  contains  the 
T  xy 

covariance  information  between  x  and  y  (x) .  If  we  let  the  vector  x  be  com¬ 
prised  of  the  i-jth  cells  (e.g.,  j  is  fixed  and  i  runs  from  I  through  I-td; 
d  is  the  dimension  of  x)  of  the  reference  array,  and.  let  y(x)  be  the  k-1 
cells,  then  A^  has  elements  given  by 

2  -/(i-k)2  +  (j-l)2/L 

T 


We  also  know  the  main  diagonal  terms  are  given  by 


2  -t/L 
dT  e  T 


Now  we  let  z  =  [x  j  ylx)]  be  the  2d-dimensional  vector  formed  by  con¬ 
catenating  x  and  y(x).  Then  the  joint  distribution  of  x  and  y(x)  is  given 
by 

,”“1/2  ^  — J  rp 

Pjj,y(5l'52)  =  P5(a)  =  ( 2it)  |Az|  exp  {-  -  5Az  5"} 


where  5^  and  a2  are  the  "halves"  of  a  corresponding  to  values  of  x  and 

y(r) .  In  the  special  case  that  x  is  large  enough  A  is  essentially  a  null 

xy 

matrix,  and  A^  has  a  block-diagonal  structure.  Then 

p_(a)  =  p-(ai)p-(a2) 
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THE  DISTRIBUTION  ON  THE  DIFFERENCE  VECTOR 

The  match  functions  defined  by  both  the  MSD  and  the  MAD  classifiers 
involve  operations  on  the  components  of  a  difference  vector.  We  now 
want  to  compute  this  difference  vector  and  determine  its  underlying 
probability  density  function. 

When  TEKCOM  samples  terrain  it  is,  of  course,  sampling  the  true  fix 
point,  x.  This  observation  is  perturbed  by  additive  independent  noise. 

The  observation  is,  therefore 

r  =  x  +  n 

Since  x  and  n  are  independent,  the  distribution  on  £  is  easily  computed 
(e.g.,  via  characteristic  functions)  to  be 

p-(a)  =  (2tt)  d//2|Arl  ^  exp  {-  j  5Ar  a1} 

where  .  A  4  A„+  A„ 

r  T  N 

The  difference  vector  that  the  classifier  uses  is  then 

d  4  r  -  y(x) 

In  the  special  case  that  x  and  y(x)  are  uncorrelated  (e.g.,  large  enough 
t)  the  distribution  on  d  is  as  easy  to  obtain  as  the  distribution  on  r. 

p-(a)  =  (2tt)  d//2|Ad|  exp  {-  j  a/i^  aT  } 

where  A,  =  A  +  A  =  2Am  +  A  . 

d  r  T  T  n 

The  more  general  case  of  correlated  x  and  y(x)  is  harder  to  analyze, 
and  is  not  considered  herein. 

THE  MSD  CLASSIFIER 

In  this  section  we  want  to  compute  the  statistics  on  the  match  function 
for  the  MSD  classifier,  and  then  compute  the  approximate  probability  of  cor¬ 
rect  assignment  according  to  equation  1,  above. 
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The  match  function  for  the  MSD  classifier  is  given  by 

d-1 

m  =  £  (r.  -  y. ) 2  =  d  •  d  =  |d|2 

i=0  1  1 

We  will  use  the  following  notation. 

p-(a)  =  (2tt) 1 2 1  ^  exp  {-  j  aZ  5T} 

This  treats  two  special  cases: 

(1)  For  T=0,  Z=A  ,  and 

n 

(2)  for  x>Lt,  I'k 2AT+An. 

To  begin,  we  know  that  d  can  be  transformed  by  an  appropriate  orthogonal 
transformation  into  a  vector,  say  u,  that  has  the  following  properties. 

(1)  The  components  of  the  vector  u  are  independently  distributed, 
zero-mean  Gaussian  random  variables  with  different  variances  (in  general) . 

(2)  If  X?  denotes  the  variance  of  u^,  we  know  that  X?  is  the  i1"*1 
eigenvalue  of  the  covariance  matrix  E. 

(3)  The  inner  product  d*d  is  preserved  under  the  transformation  so 
that  u* u  =  d*d. 

Now  we  return  to  some  statistical  ideas.  We  know  that  if  u.  is 

x 

distributed  as  N(0,X:)  then  u?  is  distributed  as  a  Gamma  distribution, 

G(a,r) ,  with  parameters  a  =  1/2X2  and  r  =  1/2.  We  know  G(a,r)  has  the 

form  a/r(r)  (ax)r  1  e  aX(x>0)  ,  hence,  Pu2(v-)  =  (2tiX?v^) -1//2exp{-v^/2X2} 

i 

The  distribution  on  u.  has  the  characteristic  function 

(1  -  j2X?m) -1/2 

Because  the  are  independently  distributed,  the  characteristic 
function  of  the  distribution  on  u'*'u  =  m  is  given  by 

Jl1  (1  -  j2X2m)-1/2 
x=0  x 

Hence  PT(nO  =  ^  f”  e^™1*  H  (1  -  j2X?w)  ”1/2du). 

J  — <o  i-0  1 
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This  is  an  exact  expression  for  the  Gauss-Markov  terrain  model  and 
the  discrete  MSD  classifier.  Unfortunately,  it  appears  there  is  no  closed 
form  expression  for  the  inverse  Fourier  transform. 

Before  finding  an  approxiniation  for  p^ (m) ,  we  note  that  in  the  case 
A?  -  A  for  all  i,  p^tm)  is  G(1/2A,  d/2).  We  further  note  that  in  the 
case  A?  =  A  for  i  =  0,1,2, ..  .1-1,  and  A?  =  0  for  all  i  such  that  I  £,  i  id-1, 
p  >m)  is  G(1/2A,  1/2). 

Now  we  want  to  evaluate  the  A^.  Specifically,  we  want  to  know  what 
are  the  eigenvalues  of  the  covariance  matrix,  E.  We  do  not  wish  to  evalu¬ 
ate  the  A^  exactly,  but  we  do  want  a  good  approximation  to  their  behavior. 

The  approximation  comes  from  noting  that  E  is  a  Toeplitz  form  (Ref  7) , 
and  hence,  as  d  +  ®,  E  assymptotically  approaches  a  circulant  matrix. 

The  eigenvalues  of  the  Toeplitz  form  and  of  its  associated  circulant 
matrix  are  distributed  identically  in  the  limit.  Further,  the  eigenvalues 
of  a  circulant  matrix  can  be  obtained  exactly  as  the  discrete  Fourier 
transform  of  the  first  row  of  the  matrix  (Ref  7 'and  5,  p  205). 

A  Toeplitz  matrix  of  order  8,  its  associated  circulant  matrix  and 
the  error  matrix  involved  in  the  approximation  are  shown  in  Fig  2.  The 
approximation  is  very  good  for  Toeplitz  matrices  whose  first  rows  contain 
many  trailing  zeroes,  or  whose  first  rows  are  essentially  constant.  This 
corresponds  in  our  model  to  «  d  and  »  d,  respectively. 

The  behavior  of  the  eigenvalues  of  a  covariance  matrix  having  elements 

exp  {— | i  -  j|/L} 

as  approximated  by  the  Toeplicz  theory  is  shown  in  Fig  3.  Also  shown  in 
Fig  3  is  a  "pulse-approximation."  The  eigenvalues  are  approximated  a 
second  time  by  a  constant  value  over  the  first  N  large  eigenvalues,  and 
by  zero  over  the  remaining  small  eigenvalues.  The  cons  cant  and  the 
number,  N,  are  chosen  to  minimize  the  mean-square-error  of  the  pulse 
approximation . 


Figure  2.  The  circulant  approximation  to  a  Toeplitz  matrix: 

(a)  a  Toeplitz  matrix  of  order  8,  (b)  the  associated  circulant 
matrix,  and  (c)  the  error  matrix. 
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Thus,  we  approximate  p(mo)  and  p^{m)  (for  r  >  Lrf)  by  Gamma  distri¬ 
butions  having  parameters  defined  in  terms  of  the  eigenvalues  of  the 
covariance  matrices  for  the  terrain  and  the  noise.  The  eigenvalues  are 
approximately  determined  by  the  Toeplitz  theory  and  the  ,fpulse"  approxi¬ 
mation.  Computer  programs  to  implement  these  approximations  are  found 
in  Appendix  A. 

We  may  note  that  this  approximation  yields  results  that  agree  with 
Schwartz  (Ref  4) .  The  eigenvalues  of  Z  are  determined  jointly  by  the 
track  length,  d,  and  by  the  correlation  lengths,  and  L^.  The  product 
of  d  and  1^(1^)  yields  the  8  of  Ref  4.  The  behavior  of  the  Gamma  func¬ 
tions  derived  from  the  approximation  agree  with  the  behavior  of  the 
approximate  distributions  given  in  Figure  2  of  Ref  4. 

As  we  shall  see,  the  Gamma  distributions  can  be  intearated  easily  to 
obtain  the  false-fix  probabilities  for  TERCOM. 

THE  MAD  CLASSIFIER 

The  match  function  for  the  MSD  classifier  is  given  by 


d-1  d-1 

m  =  l  |r±  -  yj  =  l  IdJ 


i=0 


i-0 


where  the  notation  is  the  same  as  in  previous  sections.  This  match  func¬ 
tion  is  often  referred  to  as  the  first  Minkowski  metric  because  it  cor¬ 
responds  to  K=1  in  the  following  definition  of  the  Minkowski  metrics. 


MKSiH 


a-i 

l Jri  -  yJ 


i=0 


We  may  note  that  K=2  corresponds  to  the  Euclidean  metric  found  in  the  MSD 
classifier  discussed  in  the  preceding  section. 

Because  the  Mj  metric  is  not  preserved  under  an  orthogonal  trans¬ 
formation,  it  is  very  difficult  to  obtain  the  distribution  on  m  from 
knowledge  (or  assumption)  about  the  distribution  on  d.  The  only  case 
that  is  at  all  tractable  is  the  case  of  independently  distributed  com¬ 
ponents  of  d.  This  v/ould  correspond  uo  a  diagonal  Z,  and  that  is  an 
unwarranted  assumption  for  TERCOM. 
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We  can,  however,  gain  some  insight  into  the  performance  of  the  MAD 
classifier  by  considering  a  theorem  in  pattern  recognition  that  states 
that  any  classification  procedure  is  optimal  (in  a  Bayes  sense)  for  some 
distribution  on  the  patterns.  A  meaningful  question  to  ask  about  Mj 
nearest-prototype  classifiers  is,  therefore,  for  what  distribution  of  r 
conditioned  on  y(t)  are  such  classifiers  optimal? 

Toussaint  (Ref  2)  has  shown  that  minimizing  distances  is  equiva¬ 
lent  to  maximizing  the  Laplacian  distribution.  To  make  this  clear,  let 
p(?|y(r))  be  the  conditional  probability  of  receiving  (observing)  r  when  we 
are  truly  at  location  y(t).  As  before,  ?  and  y(t)  are  d-dimensional 
ve-'tors.  If 


p(r|y(T) ) 


_1_ 

d 

a 


d-1 
II  exp 
i=0 


[- 


d-1 
II  exp{ 
i=0 


(c  a  constant) ,  then  assigning  r  to  y(i)  on  the  basis  of  nearest  Mj  dis¬ 
tance  is  an  optimal  (Bayes)  procedure.  The  extent  to  which  the  distribu¬ 
tion  of  real-world  noise  approximates  the  Laplacian  distribution  is, 
therefore,  an  indicator  of  how  close  to  optimal  the  TERCOM  classifier  is. 
While  this  model  of  the  noise  may  not  be  too  unrealistic,  it  does  not  help 
us  make  numerical  predictions  about  the  classifier's  performance.  It  does, 
however,  suggest  that  MAD  performance  and  MSD  performance  will  be  approxi¬ 
mately  the  same  when  the  noise  process  can  be  modeled  equally  well  by  the 
eaplacian  and  the  Gauss-Markov  processes.  Some  additional  points  to  con¬ 
sider  regarding  M^  classifiers  follow. 

The  Laplacian  distribution  is  not  "rotational ly"  symmetric,  so  the 
decision  boundaries  implemented  by  M^  depend  on  the  choice  of  a  coordinate 
reference  system.  This  may  be  a  distinct  disadvantage  when  modeling  a 
physical  process. 

The  distance  has  an  unusual  "instability"  that  serves  to  point  out 
both  the  coordinate  dependence  referred  to  above,  and  the  fact  that  trans¬ 
forms  that  are  usually  considered  isometries  do  not  preserve  M^  distances. 
See  Figure  4. 
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Figure  4.  Decision  boundaries  established 
by  an  MAD  classifier.  The  shaded  areas 
are  equidistant  from  yj  and  y 2  in  the 
special  case  shown. 

The  MAD  classifier  is  harder  to  analyze  mathematically  than  a  classi¬ 
fier  based  on  Euclidean  distance.  This  is  so  for  two  reasons:  (1)  An  MAD 
classifier  has  piecewise  linear  (e.g.  nonlinear)  decision  boundaries, 
whereas  the  Euclidean  classifier  has  linear  decision  boundaries.  (2)  The 
Gaussian  distribution  for  which  the  Euclidean  classifier  is  optimal  has 
been  extensively  studied  in  the  literature  due  to  the  tractability  of  its 
functional  form,  whereas  the  Laplacian  distribution  does  not  enjoy  such 
rich  development.  However,  as  Figure  5  shows,  the  difference  between  the 
two  classifiers  is  not  great  in  the  local  regions  of  the  signal  space  that 
are  "near"  the  signals  (prototypes).  Hence,  an  analysis  of  the  Euclidean 
classifier  may  produce  a  good  approximate  analysis  of  the  MAD  classifier. 
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PREDICTION  OF  FALSE-FIX  ERROR 

At  this  point  we  want  to  evaluate  Eq  1  by  using  the  Gamma  approxima¬ 
tion  discussed  earlier.  We  thereto-  •»  assume  that 


p(mc)  = 


2ozr(i) 

« 


_5HL 

2o2 

H 


1-1 


exp  {-mo/202} 


K-l 


and 


<»)  =  ^2777n  exp 


Pt>L  "  2ozr(K)  '20 


where  o2,  o2,  I,  and  K  are  the  parameters  involved  in  the  approximation. 
(See  Appendix  A  for  a  computer  program  to  find  these  once  the  terrain  and 
noise  processes  are  specified.) 

Thus ,  „ 


P(mo) 


fm  0 

p  (m)  dm  dn»o 
0  t>lt 

il-l 


-  fo  (^)  exp  (-mo/2°2»1 

■  r  woo 

By  changing  variables,  u  =  m/2o2  and  v  =  m0/2o2,  we  rewrite  this  as 
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1  i-l  -v  ,  -(a2v/o2) 

o  TUT v  e  1'e  N 


K-l  (o^v/cj2)K-1-i  ' 
i=0  <*-!-»  1 


K-l  fa2/a2)K-1-i 


i  y  iazii_  r  i-i+K-i-i  —  die 

r  (I)  >0  (K-l-i)  !  JoV 


2/o2)v 
r  dv 


K-l  (o2/a2)K-1-1 


_ 1_  p  '  # 

~  r<i)  (K-l-i)! 


(I+K-i-2) i 
, . ,  2  /  2\ I+K-i-1 

(l+oya^) 


.  .  (  j  — 1*4" 

with  R=o2/o2  and  ^  a  binomial  coefficient.  Recall  that  I  and  K  are 
half  the  number  of  nonzero  eigenvalues  of  and  of  2^T+AN>  respectively, 
as  derived  from  the  approximation. 

Appendix  A  also  gives  a  program  to  compute  this  finite  sum.  A 
recursion  formu1 i  is  used  to  avoid  an  overflow  problem  that  might  other¬ 
wise  arise  from  the  factorials  involved. 

Once  the  sum  is  evaluated,  En  1  is  easy  to  compute.  If  the  sum 
evaluates  to  a  quantity  less  than  one-half,  the  Prob  {mo  <  all  m}  is  set 
to  zero. 

We  return  to  a  discussion  of  the  predictions  from  this  model  after 
we  discuss  the  simulation.  Ac  that  time  we  can  compare  the  results  of 
the  theory  and  the  simulation. 
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TERCOM  SIMULATION 


GENERAL  CONSIDERATIONS 

For  a  theoretical  analysis  of  TERCOM  performance  to  be  useful,  one 
must  have  confidence  that  most  of  the  results  predicted  by  this  analysis 
can  be  experimentally  verified.  While  the  ultimate  test  is  an  actual 
flight  test  over  real  terrain,  a  computer  simulation  of  a  Gaussian  ter¬ 
rain  model  has  proven  very  enlightening.  Samples  of  Gaussian  terrain 
were  generated  by  a  digital  computer  program.  TERCOM  runs  across  these 
terrain  samples  were  then  simulated  for  both  MSD  and  MAD  classifiers  using 
an  uncorrelated  Gaussian  process  to  simulate  system  noise  from  all  causes. 
A  tape  of  digitized  real  terrain  amplitudes  was  obtained  which  contained 
a  wide  variety  of  terrain  types  from  flat  to  mountainous.  TERCOM  runs 
were  simulated  across  these  samples,  again  using  the  uncorrelated  Gaussian 
process  as  system  noise.  We  found  that  TERCOM  results  from  both  Gaussian 
and  real  terrain  were  almost  equivalent  when  a  certain  terrain  parameter 
ratio,  a  /'j  ,  was  the  same  for  both  types  of  terrain.  Certain  basic  con- 
elusions  about  TERCOM  performance  can  be  drawn  from  these  simulations. 

OUTLINE  OF  COMPUTER  PROGRAMS 

A  reference  array  size  of  64  x  64  cells  was  chosen  for  the  simulation 
since  this  would  represent  a  realistic  size  for  an  on-board  aircraft  com¬ 
puter  system.  The  distance  between  sample  points  (cell  size)  was  left 
unspecified  for  the  Ga.issian  data,  but  the  cell  size  used  for  the  real 
terrain  samples  was  400  ft. 

The  starting  point  for  generating  the  Gaussian  data  was  a  radially 

symmetric  autocorrelation  function  of  the  form  $>(x,y)  =  expf-Zx^+y^/L} , 
where  L  is  the  correlation  length.  Small  values  of  L  indicate  rough  ter¬ 
rain  with  low  correlation  between  points,  while  large  values  of  L  indi¬ 
cate  smoother  terrain  with  a  high  correlation  between  adjacent  points.  A 
two-dimensional  Fourier  transform  of  the  autocorrelation  function  was  com¬ 
puted  using  a  fast  Fourier  transform  algorithm.  This  procedure  generated 
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tie  power  spectrum,  P2  (10^,0^),  of  the  terrain.  Generation  of  the  actual 

terraj.il,  however,  1  quires  an  inverse  transform  of  the  complex  spectrum. 

There  are  an  infinity  of  complex  spectra  for  a  given  power  spectrum,  so 

any  complex  spectrum,  4>  (o  ,u  )  =  R  (w  ,co  )  +  jl  (to  ,to  )  ,  with  the  power 
■*  x  y  exy  mxy 

spectrum  P2  (to  ,10  )  =  $*i>  would  do,  provided  it  has  conjugate  symmetry, 
x  y 

If  $  (to  ,  to  )  =  R  (to  ,10  )  +  jl  (to  ,to  )  ,  conjugate  symmetry  implies  that 
x  y  ex  y  Jmxy 

R  (to  ,to  )  =  R  {-to  ,  -to  )  and  I  (to  ,10  )  =  -I  (-to  ,-to  ).  This  is  necessary 
e  x  y  e  x  y  mxy  mxy 

and  sufficient  to  ensure  that  the  Fourier  transform  of  will  give 

a  real  function  T(x,y)  representing  the  terrain. 

Each  real  point  of  the  spectrum,  Re(tox,»o^),  is  picked  at  random  from 

a  Gaussian  distribution  with  mean  zero  and  variance  equal  co  P(to  ,to^)/3. 

x  y 

If  Ir  (to  , to  )  I  >  P(to  , to  )  a  new  number  is  picked  until  a  value 
'exy1  x  y  r 

|Re(<ox,<Oy)  |  <_  Ptio^jtOy)  is  found.  It  follows  that  the  imaginary  part 

of  the  spectrum,  I  (to  ,10  )  ,  is  given  by  >4^  (to  ,to  )  -  R  2  (to  ,to  )  .  The 
r  mxy  *  1  x  y  e  x  y 

sign  of  I  (to  ,10  )  is  chosen  at  random  to  be  +  or  -  with  probability  1/2 
mxy 

for  each  case.  When  quadrants  1  and  2  of  the  64  x  64  spectrum  array  nave 
been  computed  in  this  way,  quadrants  3  and  4  are  generated  from  1  and  2 
using  the  conjugate  symmetry  conditions.  The  Fourier  transform  of  this 
spectrum  produces  the  Gaussian  terrain  used  in  the  simulations  of  TERCOM, 
Figure  6  shows  a  number  of  amplitude  density  plots  from  this  terrain 
superimposed  on  Gaussian  density  functions  having  the  same  mean  and  vari¬ 
ance.  The  agreement  is  quite  good.  Similarly,  good  agreement  is  shown 
in  Figure  7,  where  the  distrioution  function  for  the  terrain  data  is 
plotted  on  a  special  normal  probability  versus  amplitude  graph  paper. 
Gaussian  distribution  functions  plct  as  straight  lines  on  this  type  of 
graph.  The  steeper  the  slope  the  smaller  the  variance.  The  data  is 
fitted  quite  well  by  a  straight  line.  Figure  8  shows  two  computer  plots 
of  terrain  generated  by  the  transform  method  just  described.  In  Figure  8a 
the  correlation  length,  L,  is  4  cells  while  in  Figure  8b  it  is  8  cells. 
Increased  correlation  length  has  produced  flatter,  smoother  terrain. 
Figures  6,  7  and  8  illustrate  that  the  method  produces  terrain-like 
surfaces  with  Gaussian  amplitude  densities,  and  that  terrain  roughness 
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is  controlled  to  a  certain  extent  by  L,  the  correlation  length.  The  pro¬ 
gram  (TERC  1)  listed  in  Appendix  A  is  set  up  to  generate  and  store  on  a 
permanent  file  10  such  (different)  reference  arrays  of  dimension  64  x  64. 

The  real  terrain  tape  contained  amplitudes  from  a  large  continuous 
area  sampled  into  an  array  of  1073  x  282  cells  with  a  cell  width  of  200  ft. 
We  chose  to  use  400-ft  cells  ,  sampling  every  other  point,  and  divided  the 
area  into  eight  64  x  64  arrays  representing  a  variety  of  terrain  types. 

These  were  stored  on  a  permanent  file  in  the  same  format  as  the  Gaussian 
terrain,  so  both  could  be  used  as  input  to  the  TERCOM  simulation  program. 
Some  of  the  real  terrain  samples  are  shown  in  Figures  9  and  10.  Amplitude 
distribution  and  density  functions  were  computed  from  the  entire  1073  x  282 
array.  The  distribution  function  is  shown  in  Figure  11  plotted  on  a  normal 
probability  graph.  The  primary  deviation  from  a  straight  line  is  below  .05 
probability.  Thus  5%  of  the  amplitudes  deviate  considerably  from  Gaussian 
behavior.  The  Gaussian  density  function  and  an  amplitude  histogram  of  the 
real  terrain  cure  shown  in  Figure  12.  The  real  terrain  histogram  is  nar¬ 
rower  in  the  peak  and  wider  in  the  skirts  than  a  Gaussian  density.  We 
show,  however,  that  the  Gaussian  approximation  is  accurate  enough  for 
TERCOM  performance  estimation  purposes. 

SIMULATION  OF  TERCOM  SYSTEM  NOISE 

Simulation  of  the  TERCOM  system  required  testing  with  many  levels  of 
additive  Gaussian  noise.  The  Fourier  transform  method  used  to  generate 
the  terrain  could  also  have  been  used  to  generate  this  system  noise  if 
the  computer  memory  requirements  for  transforming  a  64  x  64  array  were 
not  so  large.  Computer  turnaround  time  on  our  CDC  6600  system  is  largely 
determined  by  memory  allocation,  so  a  more  efficient  method  of  generating 
correlated  Gaussian  noise  was  sought.  Moshman  (Ref  3)  has  demonstrated  a 
method  for  generating  a  one-dimensional  string  of  Gaussian  random  numbers 
with  any  desired  correlation  coefficient.  This  method  was  expanded  to 
two  dimensions  for  generating  a  64  x  64  array  of  correlated,  normally 
distributed  amplitudes.  The  method  is  illustrated  in  Appendix  B.  Noise 
characteristics  may  be  specified  by  correlation  coefficient  or  correlation 
length.  Both  terrain  and  noise  are  generated  from  Gaussian  random  processes 
of  mean  zero,  but  then  have  independent  variances  and  correlation  lengths. 
Most  of  the  situations  studied  used  noise  with  zero  correlation,  but  the 
option  exists  in  the  TERCOM  simulation  program  to  use  correlated  noise. 


-,<rMEAN+l<r 
RELATIVE  AMPLITUDE 


Figure  12.  Real  terrain  amplitude  histogram  compared  with  a  .normal 
density  curve  having  the  same  rean  and  variance . 


SIMULATION  PERFORMANCE  ON  GAUSSIAN  TERRAIN 

Performance  of  the  TERCOM  system  on  the  Gaussian  terrain  wi.'l  be  dis¬ 
cussed  in  some  detail.  It  will  be  seen  later  that  TERCOM  performance  on 
real  terrain  seems  to  be  somewhat  different,  but  Gaussian  terrain  per¬ 
formance  can  be  transformed  into  a  close  approximation  to  real  terrain 
performance.  Therefore,  the  conclusions  made  in  this  report  with  regard 
to  the  Gaussian  terrain  will  be  seen  to  hold  also  for  -eal  terrain. 

We  have  assumed  that  the  form  of  additive  noise  ic  Gaussian  in  the 
absence  of  any  other  evidence.  Previous  reports  by  LTV  E-Systems  have 
indicated  that  a  large  source  of  noise  is  the  error  made  in  transferring 
measurements  from  maps  or  aerial  photographs  to  the  digitized  on-board 
memory  of  the  TERCOM  system.  Such  errors  are  likely  to  be  Gaussian  in 
nature.  We  have  lumped  this  together  with  any  sensor  noise  into  one 
Gaussian  function  in  our  model.  Total  system  noise  is  varied  by  changing 

the  variance  (o2  )  of  this  noise  function. 

N 

The  TERCOM  system  in  an  aircraft  scans  a  strip  of  terrain  over  which 
it  is  flying,  recording  and  digitizing  terrain  altitudes  at  regular  inter¬ 
vals  until  some  predetermined  number  (d)  have  been  recorded.  The  on-board 
memory  contains  a  digitized  representation  of  the  surrounding  terrain  area 
over  which  the  aircraft  is  expected  to  be  flying  when  sampling  starts. 

This  digitized  terrain  is  in  the  form  of  an  M  by  N  matrix  in  which  both 
M  and  N  are  usually  larger  than  d.  We  will  assume  that  the  aircraft  is 
flying  over  the  represented  terrain  in  a  columnwise  direction.  The  dis¬ 
tance  will  then  be  computed  between  the  d-dimensional  scanned  vector  and 
all  the  column  subvectors  of  size  d  contained  in  the  array.  These 
(M-d+1) N  distances  are  computed  by  either  a  mean  square  distance  algo¬ 
rithm  (MSD)  or  a  mean  absolute  difference  algorithm  (MAD) .  In  each  case 
the  array  vector  with  minimum  distance  from  the  scanned  vector  is  picked 
as  a  match  and  the  navigation  system  uses  the  coordinates  of  this  match 
to  fix  its  position. 

Our  computer  simulation  used  a  64  x  64  array  of  amplitudes  to  repre¬ 
sent  the  terrain.  This  same  array,  with  Gaussian  noise  added,  then  repre¬ 
sented  the  on-board  memory.  Scanning  track  length  for  most  cases  was 
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chosen  to  be  48  cells.  We  assumed  that  scanning  would  start  inside  the 
specified  area  and  be  completed  within  this  area.  No  study  was  made  of 
the  situation  where  scanned  tracks  started  or  stopped  outside  the  area  of 
interest.  Thus,  in  our  case  M  and  N  were  64;  track  starting  coordinates 
were  1  to  64  in  the  N  direction  and  1  to  3  7  in  the  M  direction 
(M  -  d  +  1  =  64  -  48  +  1  =  17) .  For  each  terrain  sample,  10  such  start¬ 
ing  points  were  chosen  from  a  uniform  distribution  over  the  64  x  17  area. 

A  single  computer  run  produced  10  different  terrains  with  10  scanned  tracks 
per  terrain,  giving  100  trials  to  bhe  TERCOM  system  at  a  given  signal  to 
noise  ratio.  The  signal  to  noise  ratio  was  specified  as  where  a T  is 

the  standard  deviation  of  the  terrain  and  a  is  the  standard  deviation  of 

N 

the  noise. 

TERCOM  performance  using  both  MSD  and  MAD  classifiers  is  shown,  for 
Gaussian  terrain,  in  Figure  13.  Here  terrain  correlation  length  is  four 
cells  and  system  noise  is  uncorrelated.  The  curves  represent  the  mean  of 
400  TERCOM  trials  at  each  signal  to  noise  ratio.  One  striking  feature  of 
these  curves  is  that  system  performance  never  reaches  100%  correct  identi¬ 
fication.  The  average  miss  distance  also  remains  essentially  constant. 

This  failure  to  reach  100%  correct  at  high  signal  to  noise  ratios  follows 
from  the  method  of  "mean  removal"  used  by  TERCOM  to  normalize  the  terrain 
data. 

Mean  removal  in  the  present  TERCOM  system  as  carried  out  as  follows. 

(1)  The  mean  of  the  stored  64  x  64  array  is  computed  and  subtracted 
from  each  element  of  the  array  to  assure  that  the  stored  array  has  zero 
mean. 

(2)  The  mean  of  the  scanned  strip  is  computed  and  subtracted  from 
each  component  of  the  strip  vector  to  assure  that  the  scanned  strip  has 
mean  zero. 

(3)  The  zero  mean  strip  vector  is  compared  with  equal  length  strips 
in  the  zero  mean  stored  array  using  the  MSD  or  MAD  algorithm. 
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Figure  13.  TERCOM  performance  on  Gaussian  terrain 
with  uncorrelated  Gaussian  noise. 
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The  problem  with  this  method  is  in  step  1  above.  The  system  computes 
the  distance  between  a  d-dimensional  zero  mean  scanned  strip  vector  and  a 
d-dimensional  vector  picked  from  an  array  with  mean  zero  but  with  dimen¬ 
sionality  many  times  larger  than  d.  A  small  sample  from  an  array  of  points 
with  mean  zero  will  in  general  have  nonzero  mean.  Thus,  the  distance 
between  the  two  vectors  will  almost  never  be  zero,  even  in  the  zero  noise 
case.  This  problem  can  be  remedied  by  simply  removing  the  mean  of  each 
strip  vector  chosen  from  the  array  before  the  distance  from  the  zero  mean 
scanned  strip  is  computed.  This  alteration  gives  the  results  shown  in 
Figure  14,  where  the  terrain  used  is  the  same  as  that  which  produced 
Figure  13.  The  percentage  of  correct  identifications  rises  quickly  to 
100%,  and  the  average  miss  distance  falls  quickly  to  zero.  We  see  that 
the  simulation  has  pointed  out  a  procedural  error  not  too  obvious  from 
purely  theoretical  studies. 

A  number  of  computer  runs  were  made  to  check  out  various  effects  such 
as  correlated  noise  and  variations  in  track  length.  These  runs  were  made 
with  the  original  TERCOM  system  which  could  not  reach  100%  correct  per¬ 
formance.  The  curves  have  generally  the  same  shape  in  all  cases.  Cor¬ 
related  noise  with  correlation  length  the  same  as  the  terrain  causes  a 
noticeable  drop  in  performance,  but  the  average  miss  distance  is  reduced 
somewhat.  These  results  are  shown  in  Figure  15.  The  percentage  of  cor¬ 
rect  identifications  was  maximum  for  uncorrelated  system  noise.  TERCOM 
performance  with  this  type  of  noise,  using  both  MSD  and  MAD  classifiers, 
is  shown  in  Figure  16a.  Here  the  percentage  of  correct  identifications 
is  shown  as  a  function  of  track  length  and  terrain  correlation  length. 

In  all  cases  MSD  is  slightly  superior  to  MAD  (this  is  to  be  expected  - 
we  forced  the  process  to  be  Gaussian) .  Average  miss  distance  as  a  func¬ 
tion  of  track  length  and  terrain  correlation  length  is  shown  in  Figure  16i . 
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Figure  15.  TERCOM  performance 
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One  more  study  involving  terrain  correlation  length  was  done  using 
the  improved  TERCOM  system.  Performance  curves  were  plotted  for  increas¬ 
ing  terrain  correlation  lengths,  and  some  of  the  results  are  shown  in 
Figure  17.  We  see  that  the  curves  approach  an  asymptote  as  the  correla¬ 
tion  length  is  increased.  This  asymptotic  curve  is  important  in  interpret¬ 
ing  the  results  of  the  next  section,  where  TERCOM  performance  on  real 
terrain  is  considered. 
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Figure  17.  TERCOM  performance  on  Gaussian  terrain,  showing  asymptotic 
behavior  with  increasing  correlation  length.  Only  the  MSD  classifier 
with  improved  mean  removal  is  shown. 
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SIMULATION  PERFORMANCE  ON  REAL  TERRAIN 

Each  of  the  eight  real  terrain  patterns  was  used  for  40  TERCOM  runs 
at  each  of  seven  signal  to  noise  ratios.  The  results  of  these  runs  are 
shown  as  the  data  points  of  Figure  18.  Each  point  represents  the  per¬ 
centage  of  correct  matches  in  40  TERCOM  runs. 

Note  that  all  of  the  points  lie  below  the  solid  curve  in  the  figure. 
This  curve  represents  the  performance  asymptote  of  Figure  17  approached 
by  TERCOM  trials  on  Gaussian  terrain  as  correlation  length  is  increased. 
Thus  we  see  a  wide  variety  of  curves  for  TERCOM  applied  to  real  terrain, 
all  lying  below  the  lowest  possible  performance  curve  produced  on  Gaussian 
terrain.  Examination  of  computer  plots  of  real  terrain,  such  as  those 
shown  in  Figures  9  and  10,  revealed  a  correlation  between  initial  slope 
of  the  performance  curves  and  the  amount  of  flat  area  in  the  terrain  over 
which  this  performance  was  computed.  Pattern  1  was  completely  flat. 
Patterns  2,  3,  and  4  showed  increasing  amounts  of  mountainous  area. 
Patterns  5,  6,  7  and  8,  which  cluster  together,  were  almost  all  completely 
mountainous.  These  last  four  are  also  the  curves  nearest  the  asymptotic 
performance  curve  for  the  Gaussian  terrain. 

A  quick  analysis  of  the  TERCOM  algorithm  will  make  it  obvious  that 
TERCOM  cannot  operate  successfully  on  a  nearly  flat  or  planar  area,  and 
we  should  expect  more  mistakes  on  terrain  containing  flat  areas.  However, 
the  data  for  oT  on  the  right  of  Figure  18  show  that  there  is  no  apparent 
correlation  between  oT  and  the  amount  of  flat  or  planar  area  present  in 
the  pattern.  An  area  may  be  planar  but  have  non-zero  slope  and  conse¬ 
quently  a  fair  sized  0T.  Thus  a T  and  the  ratio  <?T/oN  are  not  the  correct 
parameters  to  characterize  performance  over  real  terrain.  Some  measure 
of  terrain  roughness  seems  to  be  needed.  Ideally,  it  seems  that  this 
should  be  some  sort  of  second  derivative  or  second  difference.  A  simpler 
measure,  which  seems  to  work  well  for  real  terrain  is 
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This  parameter  gives  a  measure  of  the  average  difference  between  adjacent 

sample  points  in  the  direction  of  the  expected  flight  path.  The  ratio 

o  /o  was  found  to  be  especially  useful  in  characterizing  the  terrain. 

A  small  value  of  a  /a  indicates  smooth  terrain  with  small  variation 
z  T 

between  sample  points,  but  possible  large  slow  fluctuations  in  amplitude 
over  tn*.  ,nole  area.  Large  0Z/0T  indicates  variation  between  adjacent 
sample  points  is  large  compared  to  the  overall  amplitude  variation  in 
the  area. 

The  ratio  °Z/°T  was  computed  for  the  Gaussian  terrain.  Because 

TERCQM  performance  on  Gaussian  terrain  approached  the  asymptote  of 

Figure  18,  the  °Z/°T  ratio  for  Gaussian  terrain  asymptotically  _ 

0.3.  This  was  higher  than  the  a  /a  values  for  real  terrain,  which  ranged 

Z  T 

from  .06  for  very  flat  areas  to  .22  for  completely  mountainous  areas.  The 

real  terrain  data  of  Figure  18  was  then  replotted  in  terms  of  a  /a  ,  where 

Z  N 

took  on  the  values  of  a^/ 1.5,  o^/ 2.5,  a^/3,  a^/ 4,  a^/7  and  o^/ 10.  This 

replotted  data  is  shown  in  Figure  19;  it  exhibits  a  very  nice  clustering 

of  the  data  points.  Three  distinct  groupings  of  curves  are  evident,  and 

there  is  correlation  between  initial  slope  and  02/crT  as  can  be  seen  on 

the  right  of  Figure  19.  Further,  all  curves  seem  to  originate  near  the 

» 

same  value  of  c  /a  . 

z  N 

Fiat  areas  were  then  introduced  into  the  Gaussian  terrain  in  the  hope 
that  this  would  decrease  their  a2/0T  below  the  asymptotic  value  of  0.3. 

The  flat  areas  were  introduced  by  simply  setting  entire  columns  of  the 
64  x  64  matrix  equal  to  zero.  This  had  the  desired  effect,  and  three  modi¬ 
fied  Gaussian  arrays  were  generated.  The  number  of  zeroed  columns  and 

corresponding  o  /a  values  were:  2  columns  zeroed,  a  /a  =  .22;  16  columns 
Z  X  Z  X 

zeroed,  »  /a  =  .10;  32  columns  zeroed,  o  /o„  =  .06.  This  modification 
z  T  z  T 

was  done  with  Gaussian  terrain  of  correlation  length  1024  cells. 
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Figure  19.  Correct  match  scores  for  MSD  classifier  over  real  terrain 
as  a  function  of  °z/aN*  The  data  show  the  emergence  of  three  groups. 

The  percentage  of  correct  scores  of  a  TERCOM  classifier  using  the  KSD 

algorithm  on  partially  flattened  Gaussian  terrain  are  shown  in  Figures  20 

and  21,  along  with  the  appropriate  real  terrain  data  points.  Note  that 

for  all  but  one  case  with  very  low  o  /a  ,  the  Gaussian  model  underesti- 

z  T 

mates  the  real  terrain  scores  by  an  average  of  about  10%.  Thus,  this 
simply  modified  Gaussian  terrain  model  gives  a  somewhat  pessimistic  pre¬ 
diction  of  real  terrain  TERCOM  performance. 
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Figure  21.  Correct  match  scores  for  MSD  classifier  over  Gaussian  and  real  terrain 
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DISTRIBUTION  OF  MISS  DISTANCES 

Typical  miss  distance  distributions  for  two  real  ..nd  two  Gaussian 
terrain  areas  are  compared  in  Figures  22  and  23.  The  values  of  ®2/oT 
for  the  terrains  are  comparable  in  each  case.  The  bin  widths  in  the 
histograms  are  5  cells,  corresponding  to  2000  ft  on  the  real  terrain. 
MSD  and  MAD  distributions  are  very  close  in  each  case  with  MSD  usually 
marginally  better  in  the  percentage  of  correct  identifications.  If 
2000- ft  errors  are  acceptable,  which  is  unlikely,  the  number  of  correct 
identifications  increase  considerably  for  most  signal-to-noise  ratios. 
From  all  the  results,  however,  completely  mountainous  terrain  with 
oz/oN  >_  1.4  gives  almost  perfect  results  with  MSD  or  MAD. 


CONCLUSIONS  AND  RECOMMENDATIONS 


COMPARISON  OF  THEORY  AND  SIMULATION 

The  results  of  the  theory  are  compared  to  the  results  of  the  simula¬ 
tion  in  Figures  24  through  27.  The  comparison  is  not  direct  because  the 
theory  predicts  how  often  the  false-fix  distance  will  be  less  than  or 
equc 1  to  the  terrain  correlation  length,  whereas  the  simulations  count 
a  decision  correct  only  when  a  correct  match  is  made.  To  make  a  direct 
comparison,  we  would  have  to  shift  the  simulation  results  upwards  or  the 
theoretical  results  downwards.  The  amount  of  this  shift  is  not  hard  to 
compute  -  either  in  theory  or  in  simulation  -•  but  it  is  a  computation 
that  we  omit  in  the  interest  of  expediency. 

In  Figure  24  we  see  two  things:  the  effect  of  increasing  track  length, 
and  the  effect  of  using  the  present  TERCOM  mean-removal  scheme.  In  view- 
Jhc  unuGi.  c  l.laI  ^>iCuxC  i»xGns ,  we  must  recall  that  the  Toeplxtr  approxi¬ 
mation  is  best  for  <<  d  or  L^  >>  d  (in  the  figures,  which  are  computer 
drawn,  N  =  d,  XT  =  L^,  and  XN  =  L^) .  For  signal-to-noise  ratios  greater 
than  about  one,  the  theory  and  the  simulation  agree  that  performance 
improves  with  increased  track  length.  In  Figure  25  the  agreement  between 
theory  and  simulation  is  quite  good.  We  again  see  the  beneficial  effect 
of  increasing  the  track  length.  A  comparison  of  Figures  24  and  25  further 
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shows  that  performance  will  improve  slightly  in  the  case  where  the  noise 
has  the  same  correlation  length  as  the  terrain.  This  same  effect  is 
observable  in  Figure  26,  where  it  is  presented  explicitly.  Again,  the 
agreement  between  theory  and  simulation  is  quite  goc*  Finally,  in 
Figure  27  we  see  that,  according  to  the  simulation,  performance  decreases 
moiiotonically  as  the  terrain  correlation  length  increases  (other  parame¬ 
ters  remaining  constant) .  The  theory  agrees ,  except  that  the  drop  seems 
to  be  too  great  for  intermediate  values  of  this  results  in  an  apparent 
increase  in  performance  as  L^,  increases  from  4  to,  say,  32.  The  reasons 
for  this  apparent  discrepancy  are  believed  to  be  due  first  to  the  Toeplitz 
approximation,  and  second  to  the  scoring  procedure.  In  this  last  state¬ 
ment  we  refer  to  the  previously  mentioned,  indirect  nature  of  the  compari¬ 
son  between  theory  and  simulation.  However,  it  seems  entirely  appropriate 
to  conclude  that  performance  decreases  with  increasing  terrain  correlation 
length. 
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DISTANCE  OF  MATCH  FROM  CORRECT  IDENTIFICATION  ( IN  CELLS) 

Figure  22.  Comparison  of  TERCOM  performance  on  Gaussian  terrain  (left  column) 
and  real  terrain  (right  column).  Open  bars  represent  the  MAD  classifier; 
solid  bars  represent  the  MSD  classifier;  one  cell  equals  400  feet  on  the 
real  terrain. 
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PERCENT  OF  TERCOM  MATCHES 


GAUSSIAN  TERRAIN  CTz/<Tr  =  . 22 


REAL  TERRAIN  <Tt/CTr=  .20 
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DISTANCE  OF  MATCH  FROM  CORRECT  IDENTIFICATION  ( IN  CELLS) 


Figure  23.  Comparison  of  TERCOM  performance  on  Gaussian  terrain  (left  column) 
and  real  terrain  (right  column) .  Open  bars  represent  the  MAD  classifier; 
solid  bars  represent  the  MSD  classifier;  on,  cell  equals  400  feet  on  the 
real  terrain. 
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AREAS  OF  FUTURE  RESEARCH 

The  results  of  this  study  partially  validate  a  test  for  TERCOM  per¬ 
formance  dependent  only  on  three  parameters:  o  and  u  (which  sunmarize 

Z  T 

the  terrain  correlation  information)  from  the  terrain  and  o„  from  the 

N 

TERCOM  system  (including  all  sources  of  noise) .  We  envision  the  develop¬ 
ment,  with  further  research,  of  a  family  of  performance  curves.  These 
would  be  probability  of  correct  match  versus  oz/a for  a  family  of  ter¬ 
rain  parameter  ratios,  °Z/°T*  Then  we  would  not  be  required  to  simulate 
a  large  number  of  individual  flights  to  determine  the  suitability  of  a 
given  terrain  area  for  use  as  a  TERCOM  navigational  checkpoint.  Compu¬ 
tation  of  the  required  parameter  ratios  would  automatically  determine  a 
point  on  the  family  of  performance  curves  and  thereby  produce  a  value  for 
the  expected  percentage  of  correct  position  identifications. 

The  parameter  ratios  a  /o„  and  a  /a.,  will  be  useful  in  analyzing 

z  T  z  N 

TERCOM  performance,  but  we  must  now  point  out  some  limitations  of  the 
present  study.  The  simulations  we re  run,  as  we  mentioned  before t  with 
a  48-cell  TERCOM  track  over  a  64  x  64  cell  terrain  grid.  This  represents 
the  so-called  "short  track  -  long  matrix"  method  which  assumes  that  the 
inertial  guidance  system  has  enough  accuracy  to  place  the  aircraft  some¬ 
where  over  the  checkpoint  area  at  the  time  the  scan  starts.  Different 
track  lengths  and  different  array  sizes  will  change  the  TERCOM  performance 
curves  in  a  manner  predicted  by  the  theory,  and  indicated  by  the  simula¬ 
tion.  Thus  the  o  Aj„  ratio  of  1.4,  which  was  quoted  as  an  ideal  operating 
z  N 

point,  would  only  hold  for  the  specific  track  and  array  size  considered 
in  the  preceding  section  of  the  report.  We  have  not  considered  the 
"long  track  -  short  matrix"  method  at  all  in  the  simulation,  although 
the  theory  includes  this  case.  Further  simulation  studies  may  be  desired 
to  determine  performance  curves  for  a  number  cf  different  track  lengths 
and  terrain  matrix  sizes.  In  particular,  model  performance  curves  for 
the  configurations  which  have  been  used  in  flight  tests  of  TERCOM  should 
be  established,  for  it  is  comparison  with  actual  flight  test  data  that 
will  determine  the  ultimate  usefulness  of  the  Gaussian  model  and  the 
parameter  ratios  and  Successful  comparison  to  flight  test 
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results  over  a  wide  variety  of  terrain  would  show  that  TERCOM  performance 
characteristics  over  any  terrain  area  could  be  predicted  by  knowing  only 
0^  and  0^  for  the  terrain  and  for  the  system  noise. 

The  mention  of  brings  us  to  another  important  limitation  of  the 
TERCOM-terrain  model  presented  in  this  report.  We  assumed  that  the  system 
noise  was  Gaussian.  If  the  real  system  noise  is  not  well  approximated  by 
such  a  process,  flight  test  performance  data  could  differ  from  model  per¬ 
formance.  If  the  process  is  non-Gaussian,  there  might  be  differences 
between  model  and  flight  test  performance  even  if  02/aT  and  are 

<samf*  for  both.  VJe  suggest,  therefore,  a  detailed  analysis  of  the  per¬ 
tinent  noise  sources  such  as  incorrect  ground  speed,  flight  path  angular 
deviation,  altimeter  errors,  etc.  This  would  produce  greater  confidence 
in  the  results  of  this  study. 

CONCLUDING  REMARKS 

We  have  shown  that  o  /o  seems  to  be  an  adequate  description  of  real 

Z  Jt 

terrain,  and  that  Gaussian  terrain  modified  by  introducing  a  certain  per¬ 
centage  of  flat  area  closely  simulates  real  terrain. 

We  have  shown  that  the  MSD  and  MAD  classifiers  produce  almost  identi¬ 
cal  results  on  both  real  and  Gaussian  terrain.  This  indicates  that  a 
theoretical  analysis  of  the  more  tractable  MSD  classifier  is  an  adequate 
theoretical  analysis  for  the  MAD  classifier. 

We  have  shown  agreement  between  the  theory  developed  herein  and  the 
results  of  the  simulation.  However,  there  is  a  strong  assumption  in  both 
the  theory  and  the  simulation  of  stationary  statistics.  This  assumption 
was  not  valid  on  the  real  terrain  investigated.  The  ratio  ^Z/«T  seems 
to  be  insensitive  to  this  nonstationarity. 

We  conclude,  on  the  basis  of  the  theory  and  the  simulation,  that 

TERCOM  -  modified  to  include  the  improved  mean-removal  scheme  -  will 

yield  acceptable  performance  on  appropriate  terrain.  What  constitutes 

"appi  opriate"  can  be  determined  either  from  o^,  L^,  L-^,  and  d  in  the 

theoretical  framework  (assuming  the  stationarity  assumption  is  valid) , 

or  from  0  /a  and  0  /o„  in  the  simulation  framework.  There  is  no  doubt 
z  T  z  N 

that  the  quantities  o  ,  0  ,  and  o„  are  the  easiest  to  obtain. 

z  T  N 
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APPENDIX  A 


COMPUTER  PROGRAMS 


1.  Program  TERCOM  results  in  the  performance  curves  shown  in 
Figures  24  through  27. 

2.  Program  TERC1  generates  Gaussian  terrain. 

3.  Program  TER  is  the  simulation  program. 


TERCOM 


XvmJckoJu  yL  ^aoo/vOaw  •  siA'U^> 
j  \ja  (s^\  i-lsA*  c ooeuvi.ew\toe  w^oAjJL^ 
'X*oj*Ai.  0^  -%2.  \, 
l.-t(*T)^o.  v^tv^  ^  LW64W)=0* 


PROGRAM  TERCGM (INPUT ,  GUTPUT,PLOT) 

DIMENSION  EIGT ( 33)  ,  EIGNI33),  A(33) 

DIMENSION  C0V(2,64),  NN(1) 

OATA  N,  XT,  XN/32,0.,  0./ 

CALL  PLOT ( 0. » -1 • ,-3) 

CALL  PLOT ( 0. ,1.^-3) 

CALL  SYMEOL  (1.  ,6.,  .15,371-PROeAEILITY  OF  CORRECT  IDENTIFICATION, 0 
*,37) 

CALL  SYM90L(1.8,5.8, .15,26HAS  A  FUNCTION  OF  S/N  RATIO, 0. ,26) 

CALL  SYMBOL (3. 5, 3. 5, 0.15, 5H  N  =  ,0.,5) 

ZN  =  N 

CALL  NUMBER (999. 0,999. 0, 0. 15 ,ZN, 0. ,1) 

CALL  SYMPOL(3.5,3.3,C.15,5HXN  =  ,G.,5) 

CALL  NUMBER(999.0,999.0,0.15,XN,0.,i) 

CALL  SYMBOL(3.5,2.1,0.15,20HXT  IS  PARAMETER  WITH, 0. ,20) 

CALL  SYM60L(3. 5,2.90,0. 15, 21HVALUE St  0,4,16  ,0.,21) 

CALL  AXIS(0.,0.,21HSIGNAL  TO  NOISE  RATIO ,-21, 8. , 0 . , 0 . , . 5) 

CALL. AXIS!  0. ,  0  ..*J.9HPF0EABILITY  CORRECT  ,19, 5 . ,  90 . ,  C.  ,  .2) 

200  FORMAT (5X,F6.2,5X,F12.6,5X,F9. 6 , 5X , F9 . 6 , 10X , 12 , 5X , F9 . 6 , 5X , 12 , 5X,F9 
1.6) 

201  FORMAT (38H1PROFA8ILITY (CORRECT)  FOR  THE  CASE  N  =  ,I3,6H  XT  =  ,F7.3 
,  10H  AND  XN  =  ,F7.3//) 

202  FORMAT (7 X , 3HSNR , 11 X , 1HF , 13 X, 4HPT AU , 8X , 1 OHP (CORRECT) , 9X, 2HNN, 1 OX, 2H 
1CN^7X,2HNT ,10X,2HCT//) 

DO  2000  NRUN=1,3 
SNR  =  0.55 
NHALF  =  N/2  +  1 
NN (1)  =  N 
Ni  =  N  +  1 
N2  =  NHALF  ♦  1 
IF(XT.EQ.0.)GO  TO  2 
DC  11  1=1, NHALF 
COV(l,I)  =  FXP(FLOATd-I) /XT) 

11  COV (2,1)  =  C. 

OC  12  I=N2,N 

COV(l.I)  =  EXPIFLOAT (I-Nl) /XT) 

12  CCV (2, I)  =  C. 

CALL  FOURT  (CCV,NN, 1,-1, 0,0) 

00  7  I  =  1  »t<  *- ALF 

7  EIGT  (I)  =  COV(i,I) 

GO  TO  14 

2  CONTINUE 

DO  13  1=1, NHALF 
1?  E 1 1 j T ( I )  =  Jt 

14  CONTINUE 

IF  (XN. £Q .0 . ) GO  TO  4 
DO  4  1=1, NHALF 

CCV(1,I)  =  EXPJFLOAT (1-I1/XN) 

8  COV ( 2, I)  =  0. 

DC  9  I  =  N 2 , N 

COV(l,I)  =  EXP  (FLOAT  (I-ND/XN) 

9  CCV (2,1)  =  0. 

CALL  FOURT (CCV, NN, 1,-1 ,0,0) 

DQ  10  1=1, NHALF 

10  EIGN(T)  =  COV(l,I) 


"b  JUJYAC,  VO*&-  X.T 
viousXoJ^e 

J  C_VacK.  O  -Aua<^TV>.  . 

7  C'YVft.  Jjxjo 

V  CAW\aSo^\t,Y\  -Wadi  A*-* 

\  4o  .  TAcxa^-j 

clclc^.  XT-O- 

CWJt?.  o  COAAaSioJ^CW  . 

r  (AV*  d,  WrVyS.  c©AAaSc^-\<N\. 

\  "MYNoAj vw  4t>  J^OMSXO.— 

o<^»^XtrVixvv3-VvArN  •  '  ° 
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T^O* 


ViC-^8  ‘v-ACyL^x: 


Gn  TO  to 
4  CONTINUE 

■00  10  I=1,KMLF 

15  EIGN(I)  =1.  <J 

15  CONTINUP  7  CrtfvryAc  P*+W-O^.C/\^\Tft3-^c/A 

C6LL  APRCX  {NHALF.NN.CS,  EIC-N)  $  ~\T  * 

dr  [NT  201,  K,  >T,  XN 


Ctx.v«_  ~  0  . 


PRINT  2  J  2 
GO  1)0)  NC  TP= 1 , 85 
VC  =  2.+SNR+SNF 
00  5  I  =  l,NHAL,r 

5  All)  =  VR* EIGT 1 1 )  *  tIGN(I) 

CALL  APRCX  {NhALF , NT, f T , A ) 

R  =  GT/CN 

cotA.J  =  PTAU(NK,NT  ,P) 

IF (ERT AU ,GT . 0 . 5) ERT AL  =  0.5 
G=l.-  2.*5RTAU 
PFINT  200,  SNR,  R,  EFTAU,  E, 
pct^RrSNR*?  .  1  FO=5.*F 

re (NCTR.EO.A.OF  .NCTR.RC.6)  GC  TO  500 
IP  (NCTR.cQ.l)  CALL  PI OT ( PR NP , ° c , 3) 
call  plctifsnr,pr,2) 

•>,/  CONTINUE 

SNR=3NR+  0. 05 
1CC0  CONTINUE 

IF (N°UN. EQ  .1)  XT  =1 . 

?(■,<’  1  XT  =  4+XT 

call  plot  e 

STCP 
ENC 


•\jaj\NpxM\  "V  xuyi&L  r^\^x>6 

3  <^-<rT/5Vi 

CTeftv^ute  ^  (.TA  4.  fAji  ■ 

N N ,  CM,  NT,  CT 


^  SUVL  ^  * 

x  eiko^-Q  PO-NO^oS-aN  XT  ^ 

^  O^s  V\  XT  , 

a  ' ^tSrttvv  Wx»^  W-^ 


APPQX 

SUBROUTINE  APROX(NHALF,<,C.A) 
OI^PNSION  A (NH ALp) 

C  =  Ad) 

£n  A  x  =  All) 

K  =  1 

DC  T  1=2  , NHALF 
SUP  =  0. 
no  1  J=i,I 

1  SUP  =  SUP  +  A ( J ) 

E  =  SUB* *2 /FLOAT (I) 

IF  (E.LT.EMAX)G0  T02 
EMAX  =  E 
X  =  I 

C  =  SUM/PLCAT (I) 

?  CONTINUE 
3  CONTINUE 
RF  TURN 
EN0 


FT  AU 

FUNCTION  P T  AU (T  N  ,NT , F ) 
PI  =  P+1. 

C  =  (P/Pl) **NN 
SUM  =  1. 

IF  (NT.EO.DGO  TO  2 
A  =  i 

NT1  =  NT-1 
XNN  =  NN-1 
00  1  1  =  1,  NT1 
XI  =  I 

A  =  A* ( XNN+XI) / ( XI *R1) 

1  SUM  =  SUP  +  A 

2  PTAU  =  1.  -  G* SUP 
RETURN 

ENO 
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TERC1 


\  A. 


PROGRAM  TERCH  INPUT ,  OUTPUT  ,  TAPED 
DIMENSION  A  (2, 64,64)  ,NN (2) 

DATA  NN ,  TC , P , B/64, 64 y 4 •  ,  5 •  ,  1  •  / 
CALL  RANSET(P) 

DO  1 0 0 G  NCTR=1,10 
DO  10  1=1, 2  Vl  / y 

DO  10  J= 1, 64  '  A' 

00  10  K= 1, 64  * 

A ( I , J, K) =0. 

10  CONTINUE 

00  20  J= 1, 33 
DO  20  K=l,33 
CJ1=FL0AT  ( J) 

CK1=FL0AT(K) 


B 


C=SQRT < (CJl-i.) **2  +  (CKl-l.)  **2) 
A  ( 1 ,  J , K) =B*EXP  C-C/TC ) 

20  CONTINUE 

00  30  K= 1, 33 
DO  30  J=  2, 32 
J 1  =66- J 

A(1,J1,K)=A(1,J,K) 

30  CONTINUE 


i  'A'. 

F ',  *«  ■ 


DO  40  J= 1, 64 
DO  40  K=2, 32 
Kl=66-K 

A(1,J,K1)=AU,  J,K) 

40  CONTINUE 

CALL  FOURT  ( A ,NN , 2, , 0 , 0 )  —  Cc,  -  y 

00  60  J= 1, 64  -  ,A 

00  60  K= 1, 33 
F2=A(1,J,K) 

CALL  RAN00M(F2,RE,XI) 

A ( 1, J, K) =RE 
A  (2, J, K) =XI 
60  CONTINUE 

DO  70  J=  34 , 64 
Jl=66- J 

A ( 1 , J, 1) =A (1 ,  J1 , 1) 

A(2,J, 1) =-A (2, J1 , 1 ) 
A(1,J,33)=A(1,J1,33) 

A( 2, J. 33 ) =-A (2, Jl, 33) 

70  CONTINUE 

DO  80  K=34 ,64 
Kl=66-K 

A  ( 1, 1,  K)  =A  ( 1 , 1 ,  Kl) 

A  ( 2, 1 ,  K)  =*•  A  ( 2, 1,  Kl ) 
A(1,33,K)=A(1,33,K1) 

A { 2, 33 ,K ) =-A (2 , 33, Kl ) 

80  CONTINUE 

00  85  J=  2, 32 
DO  85  K=  2, 32 
Jl=66- J 
Kl=66-K 

A { 1, Jl ,K1) =A  ( 1 , J  ,K ) 

A ( 2, Jl , Kl) =- A ( 2, J, K) 


■  1  t-  i  V 

•:  -  rvA/ 


5 


v  tE'i'i 
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TERC1 


85  CONTINUE 

DO  90  J=2,  32 
00  90  K=2,32 
Kl=66-K 
Jl=66- J 

A(l, J,K1)=AC1, Ji,K) 

A(2,J,K1)=-A(2,J1,K) 

90  CONTINUE 

CALL  FOURT (A ,NN, 2, 1, 0 , 0) 

CALL  RANGET(P) 

100U  WRITE(l)  NCTR,  TC,  P,  B,  C(A(1,I,J),  1=1,64),  J=i,64) 
PRINT  102,  NCTR,  P 

102  FORMAT (10X,*NCTR=*I2,5X,*P=*E15.7) 

REMIND  1 

STOP 

ENO 


RANDOM 

SUBROUTINE  RANDOM ( r2 , PE, XI ) 
IF  (F2. GT .0 . ) GO  TO  1 
RE  =  0. 

XI  =  0. 

GO  TO  3 

1  F  =  SQRT  (F  2) 

ZMU=0. 

SIGMA=F/3. 

2  X=RANF (DUM ) 

Y -RANF (DUM ) 

CC=SQRT(-2.*ALOG(X) >*SIGMA 
RE=CC*COS(6.28313*Y> +ZMU 
IF  (ABS  (RE)  . GT . F ) GO  TO  2 
XI=SQRT (F2-RF**2) 

IF ( X . G  T • .5)XI=-XI 

3  RETURN 
END 


60 


TER 


6 


PROGRAM  TER  ( INPUT , OUTPUT , TAPE  1) 

OIMENSION  A (64,64) ,S(64) ,  AS  (64, 64) 

0 1 MENS  ION  RN  <  64 )  .  „  „  „  _ 

ro=o •  _  k>os~  CaR *C. LAnow  Co£fr»(i8vT 


£=.7642672 

NL=4B  - —  S £T  TRACif  LENGTH 

NBL=64 - S£r  Aff/?*y  SI2C 

NL5=NBL-NL 
BL5=PL0AT { NL5) 

CALL  RANSET(E) 

PRINT  107, E 
PRINT  110 


ANUH=SNUM=0H<50  =  0MA0=C . 

00  70  NC=1 , 1 0 

RE  AO  ( 1 )  NCT  R ,  TC ,  P  ,  8  ,  A  —  READ  75?<>At.V  PAlTE Q.'-J 
00  5  1=1,64 
00  5  J=1 ,2 
A ( I,J» =0  . 

CONTINUE 
PRINT  104, NC 
SUM10=0. 

SUM11=0. 

RA  =  2.  -  SET  ?lG.VAL  TV  VOlSf.  ?AHC 

FMEAN=0. 


IK=JK=64 


BN=FLOAT (IK*JK) 

SUM=  3 « 

8N2=FL0Af (IK-1) 

00  10  1=1, IK 

oo  n  j= i , jk 

SUM=SUH+A( I, J) 

10  CONTINUE 
AVE=SUM/BN 
SUH=). 

00  23  J= 1 , JK 
00  20  1=1, IK 
SUM=SUM* ( A ( I , J) - AVE) **  2 
20  CONTINUE 

SIGT  =  SQRT(SU«/BN)  -<?i 

SIGNO  =  SIGT/RA 
10  =  1 

IF  (RO.NE.il.)GO  TO  22 
RO 1  =  0 . 

GO  TO  25 

? 2  IF (RO. NE .1 .) GO  TO  23 
RO 1=1 . 


CovzjrA*  ,:v 

'  •  V  0  . '  A  •  ,  V 


GO  TO  25  .  .  ^vcii7(.jw  LEoGTrf 

23  SCN=-1./AL0G(R0)  -  "  •  -rouWc't JU£-Vi8<A_  . 

R01=EXP(-S0PT(2.)/SCK>  -c  .A-  •>/  Bf  rwC.£  V 
25  ACA=R0/ ( 1.  +R01)  'J'4 

AC9  =  S0RT (1  .-2.M (RO) 2) / ( R01 1 . ) ) 

CALL  RAND2  (PN1 , RN2  ,  SIG NC , FME  AN  ,  I D)~" ’ 

RN3=RN2 

DO  3  3  1=1, 6**, 2 

CALL  PANC2 (RN1,RN2,SIGNC,FMEAN,I0) 


01 


TFR 


70 


31 


52 


53 


40 

41 


51 


55 


RN(I)=R0*RN3*SGRT(i.-R0**2)*RNl 
RN ( I+l ) =RO*RN ( I ) +SGRT(1.-R0**2)*RN2 
SUM19=RN (I)**2+RN(I+i> **2+SUM10 
SUM11=RN(I)*RN(I+1)+SUM11 
AB (I»1)=A<I,1) 

A8(H-l,i>  =  A(I*i,l) +RMI+1) 

RN3=RN ( I  +1  > 

CONTINUE 
DO  31  J=2*  64 

CALL  RAND2 (PN1 , RN2 , SIGNO , FMEAN , ID) 
RN(i) =R0*RN<1) +SQRT(1 .-RO**2> * RN1 
RN (2)=ACA* (RN<1) +RN(?> ) +ACB*RN2 
SUMn  =  RN(l)**2+RN(2)  **2+SUM10 
SUMU  =  RN ( 1 ) +SN { 2 )  +  SUM1 1 
ABC1,J)=A(1,J) +RN ( i) 

A8 (2»J)=A( 2 ,  J ) +  R  N  (  2 ) 

DO  31  1=3,64,2 

CALL  RAN02 (RN 1 , RN2 , SIGNO , FMEAN , 10) 
RN  (I)  =  ACAMRN<  I'l) +  RN ( I ) ) *AC8*RN1 
RN  (H-1)  =  ACA*  <RN(  I)  +RN  (1*1)  )+ACB*RN2 
SUM10=RN (I)**2+RN(I+1) **2+SUM10 
SUM11=RN(I) +RN(Itl)+SUHll 
AB (I, J)=A( I, J) +RN( I) 
AB(I+l,J)=A(I+i, J)  +RMI  +  1) 

CONTINUE 


Nb\st 

PHIEPN 


lro •  -jw 


SUM=0. 


DO  52  J=  1 »  JK 
DO  52  1=1, IK 
SUM=SUM+  AS  (I,J) 

CONTINUE 
SUM= SUM/4096. 

00  53  J= 1 , JK 
00  53  1=1, IK 
AB(I,J)  =  AB(I,J)-SUM 
CONTINUE 

SOEVN=SORT (SUM 10/4 096. -(SUM11/ 4096. ) **?> 
IF (SOEVN.NE.C. )G0  TO  40 
SNR  =  11 . 1 1 

GO  TO  41  , 

SNR=SIGT/SOEVN  -  i  ^CT N 

PR  INTI 13  ,SDEVN,SIGT,SNR 
10  =  2 

DO  ? 0  L= 1, 10 

CALL  RAND2  (RN1 , RN? , S IGNC , FME AN , I D) 
KI1=IFIX  (RL5*RN1)  +  1 
KI2  =  KU  +  NL-1 

KJ1  =  IFIX(63.*RN2M  1  T  / 

DC  51  1=  KT  1, K I?  j  -a-  '  , C.  C}O»0' 

M= I-KI 1+  1 
S(M)=A(I,KJ1) 

CONTINUE  - 

SU«=0. 

00  55  1=1, NL 
SUN=SUM*3( I) 

CONTINUE 


62 


SUM=SUM/NL 
00  56  1=1, NL 
S(I)=S(I)-SUM 
56  CONTINUE 
NL2=64-NL+1 
SUMR=320  00  . 

SUMR2=1.E7 
DO  62  J= 1 , JK 
DC  62  K= 1 »  NL? 

SUM=T. 

SUM2  =  0  . 
su«5=o. 

00  58  1=1, NL 

II  =  I+K-i  IASD  •'AA'i 

SUM5  =  fl8(  II,  J)  +SUM5  -- 

58  CONTINUE  '  '*  ' 

SUH5=SUM5/FL0AT(NL) 

DO  60  1=1, NL 
I 1= I +K-1 

SUM=SUM+ABS(S(I)-flB(II,J)+SUH5) 

SUM2=SUM2+ (S(I)-AB(II, J) +SUM5) **2 

60  CONTINUE 

IF (SUH.6E. SUHR ) GO  TO  61 

SUMR=SUM 

KIR=K 

KJR=J 

61  IF(SUM2.GE.SUMR2)GO  70  62 
SUMR2=SUK? 

KIR2=K 
K  J,"T2  =  J 

62  CONTINUE 

OMAO=SQRT(  (FLOAT  (K  1 1-K  IR ) )  **Z  MFLO AT  (KJ1-KJR)  )**2)  +DMAD 
OMSO=SQRT( (FLOAT (KI1-KIR2) ) *  *  2  +  ( FL OAT ( KJ1-KJR2) )**2)+0MS0 
IF(KI1.E0.KIP.AN0.KJ1.E0.KJR)G0  TO  63 


IF (KI1.EG.KIR2.ANO.KJ1.EO. KJR2 ) GO  TO  66  1 

GO  TO  67 

63  IF(KI1.EO.KIR2.ANO.Kj1.EO.KJR?)GO  TO  65 
PRINT  J  C 5 , KIR , KJR, KI 1 , KJ 1 , KIR2 , K JP2 
SNUM=SNUM  <■  1  • 

GO  TO  70 

65  PRINT  106, KIR, KJR, KIP2  ,KJR2 
GO  TO  70 

h6  PRINT  108, KT 1 , KJt , KIP, KJR, KT 1 . KJ1 
ANUM-=  A  NUM*  1  . 

GO  TO  70 

67  PRINT  10<3,KI1,KJ1,KIP,KJC,KI1,KJ1,KIR?,KJR2 
ANUM  =  ANUP«-1. 

SNUM=SNUR+1. 

n  CONTINUE 

OMAO=QHAC/ ASUM 
OHSO=OMSO/SNUH 
REWIND  1 
PRINT  103, RA 
PRINT  111, CPAO,OMSO 
PRINT  112, RO 


t-  I;  -A  , 


PRINT  114, NL 
CALL  RAN GET (00) 

print  107,00 

STOP 

103  FORMAT (IX, ’SIGMA  T  /  SIGMA  N  =’,F4.1) 

104  FORMAT (IX, ’NCTR=’, 13) 

105  FORMAT  ( IX, ’CORRECT  I  .D  .  AT’  ,2  T  5 , 1  OX  ,  2  I  ' ,  *1X,  ’M I  ST  AKEN  FOR’,  213) 
l(i  6  FORMAT  (IX, ’CORRECT  1 .0  .  AT  ’  ,  2  IT  ,  1  OX  ,  ’  CORRECT  I. 0. AT’, 213) 

107  FORMAT  (IX,  ’RANDOM  NUf-  B  ER  =  ’  ,  E  15 . 7 ) 

10A  FORMAT (LX, 213, IX, ’MISTAKEN  FCR’ , 213, EX , ’COPRFCT  I. 0. AT’, 213) 

10R  FORMAT (IX, 21 3, IX,  ’MISTAKEN  FOR’ , 21 3 , EX , 2 13 , 1 X ,  ’ MIST AKFN  FOR*, 213) 

110  FORMAT (7X, ’MAO’ , 32X,’MS0’) 

111  FORMAT (IX, ’AVE.  MAO  F R ROR= ’ , 1 X , F5. 1 , 1 1 X , ’A VE .  MSO  ERROR=’ , IX , F5. 1 ) 

112  FORMAT (IX, ’R0=’,F5.?) 

113  FORMAT (60X , ’NOISE  S I C M A= ’ , E15 . 7 , IX , ’ SIGN AL  SIGM A=’ , Et 5 . 7 , 

1*S/N  =  ’,F«3. 2) 

114  FORMAT (IX, ’TRACK  LENGTH  =’,I3> 

END 


X  =  R«F  (DUM|RaN02t,,N1’RNJ’SIGN0’F,<EiN’,DI 
Y=RANF (OUM) 

IF(I0.FQ.i,G0  TO  5  -  IDM  TOfi  ^AOD.A.v  0 HTfiiGJ'ji, 

RN1  =  X  f'j*  D  •  ,v 

RN2  =  Y 
GO  TO  10 

5  CC=S0RT(-2.’AL0G(X) ) ’SIGNO 
RN1  =  CC’CCS  (6.28313’Y) fFMEAN 
RN2  =  CC’SIN  (6.28313’Y) +FMEAN 
10  RETURN 
ENO 


APPENDIX  B 
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This  can  be  extended  to  two  dimensions  by  making  S .  .  depend  on  S .  .  . • 

i,3  i~l,3 

S.  ,  .  ,  and  S.  .  , 


Si-1, j-1 


Si, j-1 


S.  ,  . 

o  1-1,3 


-OS.  . 
i,3 


we  have  then 

E(R  )  =  o2 
n 

E(S.  .S.  .)  =  O2 
i,j  1,3 

E (S .  .S.  ,  .)  =  E (S .  .S.  .  ,)  =  po2 

1,3  i“l,j  1,3  l,j~l 


E (S .  .S. 


i, j  i-1, j-1 


)  =  P ,o" 


Let  S.  .  =  w,S.  ,  .  +  w.S 


+  w_R 


i,j  1  i-l,j  2  i, j-1  "3  n 


where  again  0  <  ^  £  1,  0  £  w2  1,  0  £  w^  £  1 

and  Var  (S.  .)  =  Var  R 
i,3  n 

Var  S. .  =  E(w,S.  .  .  +  w_S.  .  ,  +  w,R  }2 
i]  1  i-l,3  2  i»3“l  3  n 

o2  =  w^  2o2  +  w 52o2  +  w^2o2  +  2w^w2p^02 


or  1  =  w^2  +  w22  +  w32  +  2w1w2p ^ 


Multiply  Eq  B2  by  S.  .  .  and  take  the  expectation  of  both  sides. 

i-i,  3 

This  gives: 


po2  =  w^o2  +  w2p1a2  +  0 


p  =  W1  +  w^ 


Multiply  Eq  B2  by  3^  ^  ^  and  take  the  expectation  of  both  sides  to  get 


po2  =  w  p  +  w2a2  +  0 


r\  ~  T.t  ^  X  u 

-  -n  -2 


Now  multiply  Eq  B4  by  p^  and  subtract  Eq  B5  from  the  result  to  get 

ppi  "  pi”i  "  Vi! 

~  P  *  -pJl“l  -  “2 
ptp^l)  =  w2(p12-l) 


p(Pl-l)  P 

(Pl^-l)  ‘  Pi  +  1 
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I 


1  1 


Similarly  we  find  that  wj  =  W£  = 


P 

Pl+1 


Mow  substituting  for  w^  and  wj  in  Eq  B3,  we  get 


1  = 


2p2 

(Pl+D2 


+ 


2PlP2 


or 


w3 


(2p2  +  2pip2) 

(Pi  +  1)2 


2P2 

(Pi  +  1) 


Thus 


S.  . 


Pi  +  1  i-l/j 


Pi 


1  Si/j-l 


+ 


2P2 

Pl  +  1 


R 

n 


When  p  =  0,  S.  .  =  R  and  wh jn  p  =  pi  =  1,  S.  .  =  -*  S.  .  ,  +  r  S.  .  , . 

1,3  n  1  13  2  1-1,3  2  1,3-1 

There  are  two  ways  to  specify  the  characteristics  of  the  correlated 

noise  we  wish  to  generate.  The  noise  may  be  characterized  by  a  correlation 

length  t  or  a  correlation  coefficient  p.  In  both  cases  an  exponential 

autocorrelation  function  is  assumed  for  the  noise. 

If  the  S.  .  are  treated  as  samples  from  a  continuous  function,  we 
1,3 

can  determine  p  and  p*  when  the  correlation  length  x  is  given.  If  the 

correlation  length  is  x  cells,  p  is  the  value  of  the  correlation  function 

one  cell  away  from  the  origin.  Thus  we  have 

_  1 
p  =  e  x 


Since  S.  ,,  S.  ,  is  1.414  cells  distant  from  S.  ., 
l-l  3-1  i/3 

_  1.414 

Pi  =  e  x 

If  the  noise  is  specified  in  terms  of  a  correlation  coefficient  p, 
we  can  determine  x  from  the  relation  x  =  l/-lnp.  Then,  Pi  can  be 
determined  as  shown  above. 


jt 

Sr 

f 
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